Beyond the Gene List

How to manipulate a gene list

The golden rule before we start Always keep a backup of the output from the pipeline, never work on the original data!

There are various tasks you might wish to perform on a gene list

• Search for our favourite gene
• Sort / Rank according to statistic or p-value
• Filter to obtain genes with particular cut-off and fold-change
• Re-order / re-arrange columns

Our recommended tools to perform such operations would be a programming language such as R / Python so that the operations can be scripted and automated. As a compromise, our Galaxy server provides various operations through an intuitive interface and allows the user to build a workflow to a chain of simple operations.

There is a special reason why you shouldn’t trust Excel as your many means of interrogating the results. Sometimes it can be unhelpful and convert gene identifers into dates. No, really!

We will only scratch the surface of what is possible with Galaxy. If you are new to the software, you can check out our course:- http://galaxycam.github.io/

Different online tools or GUIs will take different types of input, but will probably involve some combination of the operations mentioned above. Therefore, we will introduce some tools in Galaxy that should allow you to manipulate your data into the required form. Later, we will give examples of workflows for particular gene set enrichment / pathways analysis tools.

In the first generic example, we will show some common manipulations on a gene list using Galaxy:-

Common tasks in Galaxy

Upload the gene list in csv form

The gene list .csv file can be uploaded into Galaxy. In this example we want to choose the file with no p-value cut-offs applied

Convert to tabular

Before we can go further, Galaxy needs to convert the data you just uploaded into tabular form by replacing commas with a tab. This can be done using the Text Manipulation -> Convert delimiters to TAB option. Make sure that Commas is selected from the drop-down.

Filter

We can apply filtering using the menu option Filter and Sort -> Filter data on any column using simple expressions. Here we use the condition c7 <0.01 and c3 > 1.5 to require that Column 7 (adjusted p-value) is less than 0.05 and Column 3 (log fold-change) is greater than 1.5.

You could also use this menu to filter on genes on a particular chromosome or within a certain range of start and end positions

Selecting columns

If you find the number of columns in the file a bit unwieldy, you can also use the Text Manipulation -> Cut columns from a table to remove some columns, or re-order the existing ones:-

Sorting by a column

We can also sort the table by values in a particular column. e.g. the log\(_2\) fold change which is found in column number 3.

• can also sort by multiple columns
  • e.g. Chromosome, and then start

R workflow

For those keen on R, equivalent operations can be performed with the dplyr package. See our intermediate R course for details.

library(dplyr)
deTable <- read.csv("t47d_Treatment_DEA_Prog-vs-Control_all.csv")
filteredTable <- filter(deTable, padj < 0.05, log2FoldChange > 1.5)
sortedTable <- arrange(filteredTable, log2FoldChange)
reducedTable <- select(filteredTable, X,baseMean,log2FoldChange)

Gene-Ontologies and Pathways

In the early days of microarray analysis, people were happy if they got a handful of differentially-expressed genes that they could validate or follow-up. However, with later technologies (and depending on the experimental setup) we might have thousands of statistically-significant results, which no-one has the time to follow-up. Also, we might be interested in pathways / mechanisms that are altered and not just individual genes.

In this section we move towards discovering if our results are biologically significant. Are the genes that we have picked statistical flukes, or are there some commonalities.

There are two different approaches one might use, and we will highlight the theory behind both

Over-representation analysis

• “Threshold-based” methods require defintion of a statistical threshold to define list of genes to test (e.g. FDR < 0.01)

Let’s take the set of genes H2AFY, SPRY2, CIB1, WNK1, SPRED1, SPRED2. With the exception of one, they all look to have significant p-value. But is this evidence enough that the whole pathway has been disturbed?

The question we are asking here is;

“Are the number of differentially expressed (DE) genes associated with Theme X significantly greater than what we might expect by chance alone?”

We can answer this question by knowing

• the total number of DE genes
• the number of genes in your gene set
• the number of gnes in your gene set that are found to be DE
• the total number of tested genes

The formula for Fishers exact test is;

\[ p = \frac{\binom{a + b}{a}\binom{c +d}{c}}{\binom{n}{a +c}} = \frac{(a+b)!(c+d)!(a+c)!(b+d)!}{a!b!c!d!n!} \]

with :-

	In DE List	Not In DE List	Total
In Gene Set	a	b	a + b
Not in Gene Set	c	d	c + d
Total	a + c	b +d	a + b + c + d (=n)

In this first test, our genes will be grouped together according to their Gene Ontology (GO) terms:- http://www.geneontology.org/

Software for conducting a over-representation test (goseq)

We will be using goseq, which is a software package available through Bioconductor. However, rather having to write R code, we will be using the package through our institute’s Galaxy server.

This package has been specifically-developed for use with RNA-seq data. Plenty of methods have been applied to microarray data, but the assumptions might not hold for RNA-seq. For instance, goseq will adjust for the length of a gene.

Preparing the data for an over-representation test

We can start with the list that contains results for all genes; Let’s take t47d_Treatment_DEA_Prog-vs-Control_all.csv

The goal is to obtain a table with two columns; the first containing the gene identifier, and the second being TRUE or FALSE indicating whether the gene is differentially-expressed at a given cut-off (e.g. 0.05). We don’t have to be too rigourous when deciding the p-value cut-off, as we’re using it here as a means to filter the data not the end-point of our analysis (I’ve seen some people use a threshold as high as 0.1).

If you wish to see the R code (using dplyr) for manipulating our data to provide an input file for the gene set analysis, click the ‘Code’ button on the right.

library(dplyr)
deTable <- read.csv("t47d_Treatment_DEA_Prog-vs-Control_all.csv")
mutate(deTable, DE = padj < 0.05) %>% 
  mutate(DE = ifelse(is.na(DE),FALSE,DE)) %>% 
  dplyr:::select(X, DE) -> filteredTable
  write.table(filteredTable,"de-table-for-goseq.txt",quote=FALSE,row.names=FALSE)
filteredTable

Manipulating the file can also be done in Galaxy.

Galaxy Workflow

Import the csv file

Get Data -> Upload Data

Convert to tabular form

Text Manipulation -> Convert delimiters to TAB

Select Commas from drop-down menu

Remove header line

Text Manipulation -> Select last lines from a dataset

Make note of how many lines in file, N and remove N-1 last lines with this tool

Retain columns of interest

Put c1,c2,c3,c4,c5,c6,c7,c12 to retain columns 1,2,3,4,5,6,7 and 12

Text Manipulation -> Cut columns from a table

The cleaned dataset

The output from the previous step will be referred to as our cleaned dataset and we will use it later on for another type of downstream analysis.

We now need to compute if each gene is DE at 0.05 significance level

Adding an extra column

Text Manipulation -> Compute an expression on every row

Use the expression c7<0.05 to test if the adjusted p-value is less than 0.05

Which should give something like this

Extracting the columns needed for goseq

Text Manipulation -> Cut columns from a table Choose c1,c9 as the columns to cut; the ensembl ID and the column created in the previous step

Which should give the following

Running goseq

NGS:RNA-seq -> goseq

Select Genome version hg19 in this case. You could also specify which type of identifiers can be found in the first column. Ensembl ID is the default and correct for our dataset.

After a while, you should notification that your analyis has completed and you can download the results by clicking the floppy disk icon.

Things to try:

• Filter the output from the goseq analysis to gene sets with p-value less than 0.05

The differentially-expressed genes that we identified could either by up- or down-regulated. Sometimes we might want to look for over-represented genes in the up- or down-regulated genes only

• Let’s suppose we want to analyse the up-regulated genes only
• Consider what changes you would need to make to the workflow and perform the analysis

Threshold-free analysis

For these tests, we don’t have to specify a statistical threshold and use the test statistics from all genes as the input to the test. The popular Gene Set Enrichment Analysis (GSEA) method uses this approach. These tests can be used to detect differential expression for a group of genes, even when the effects are too small or there is too little data to detect the genes individually.

Let’s say that we have 10 genes “of interest”. We can look-up their statistics from the table

Attaching package: ‘limma’

The following object is masked from ‘package:BiocGenerics’:

    plotMA

Some of the statistics look high, but how do they compare to the rest of the genes? Are they towards the top of the list, or scattered randomly? This is what GSEA tries to determine. We are not making any claims about any particular gene in the list being DE, only that the set as a whole is altered.

We typically visualise the results of the analysis as a barcode plot. Firstly, all the test statistics are ranked from highest to lowest. We then draw a vertical line at every point where a gene in our list of interest is observed. We then determine where the peak of the distribution is, and consequently whether the genes are towards to top or bottom of the ranked list.

Here is the plot for our set:-

The case where we have a down-regulated pathway might look something like:-

Or in the case of no enrichment, there is no particular pattern to where the genes occur in the ranked list

The Broad institute provides a version of GSEA that can be run via a java application. The input to GSEA is matrix of expression values / counts and a file describing which covariate each sample belongs to. However, whether or not the method can be applied directly to RNA-seq data is under debate. According to the GSEA faq

The GSEA team has yet to determine whether any of these ranking statistics, originally selected for their effectiveness when used with expression data derived from DNA Microarray experiments, are appropriate for use with expression data derived from RNA-seq experiments. We hopefully will be able to devote some time to investigating this, but in the mean time, we are recommending use of the GSEAPreranked tool for conducting gene set enrichment analysis of data derived from RNA-seq experiments.

The GSEAPreranked tool mentioned is available through the Broad insitute’s GenePattern analysis suite. However, we don’t want to introduce yet another analysis suite, so we will focus on how to perform this analysis in Galaxy.

GSEA uses a set of pre-defined gene sets in their analysis which come under the following categories:-

• H hallmark gene sets
• C1 positional gene sets
• C2 curated gene sets
• C3 motif gene sets
• C4 computational gene sets
• C5 GO gene sets
• C6 oncogenic signatures
• C7 immunologic signatures

You can find out more about these gene sets on the msigdb website. However, it is important that these gene sets have Human gene symbols only.

Preparing the gene set for a GSEA analysis

For those interested, here is the R code to produce a file suitable for GSEA analysis:-

library(dplyr)
deTable <- read.csv("t47d_Treatment_DEA_Prog-vs-Control_all.csv")
deTable %>% arrange(desc(stat)) %>% 
  dplyr:::select(symbol, stat) %>% 
  filter(!is.na(symbol)) -> orderedTable
orderedTable

Now lets walkthrough how to do that in Galaxy:-

Ranking the genes

From the cleaned table created above, sort on column 5 (the test statistic)

Filter and Sort -> Sort data in ascending or descending order

Giving the output:-

We now have all the information we need, but the GSEA tool expects only two columns; gene symbol and ranked statistic

Cutting the columns

Extract columns 8 (gene symbol) followed by test statistic; c8,c5

Text Manipulation -> Cut columns from a table

Running a GSEA analysis

Using the fGSEA tool

fGSEA is a Bioconductor package that implements the pre-ranked GSEA analysis. It also happens to be a faster implementation, as described in their pre-print

You can find the tool under NGS: RNA-seq Analysis

With this tool, three different outputs are produced

• a results table
• a summary graphic of the top pathways
• enrichment plots for the top pathways

Things to try:

• Can you identify pathways that are significant at 0.05 level?
• What pathways are over- and under-represented?